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ABSTRACT. The refolding from stretched initial conformations of ubiquitin 
(PDB ID: lubq) under the quenched force is studied using the Cq-Go model and 
the Langevin dynamics. It is shown that the refolding decouples the collapse and 
folding kinetics. The force quench refolding times scale as tf ~ exp{fqAxF/kBT), 
where /y is the quench force and Axp ~ 0.96 nm is the location of the average 
transition state along the reaction coordinate given by the end-to-end distance. 
This value is close to Axp ~ 0.8 nm obtained from the force-clamp experiments 
[J. M. Fernandez and H. Li, Science 303, 1674-1678 (2004)]. The mechanical and 
thermal unfolding pathways are studied and compared with the experimental 
and all-atom simulation results in detail. The sequencing of thermal unfolding 
was found to be markedly different from the mechanical one. It is found that 
fixing the N-terminus of ubiquitin changes its mechanical unfolding pathways 
much more drastically compared to the case when the C-end is anchored. We 
obtained the distance between the native state and the transition state Axup ~ 
0.24 nm which is in reasonable agreement with the experimental data. 
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Introduction 

Deciphering the folding and unfolding pathways and free energy landscape of biomolecules 
remains a challenge in molecular biology. Traditionally, folding and unfolding are monitored 
by changing temperature or concentration of chemical denaturants. In these experiments, 
due to thermal fluctuations of initial unfolded conformations it is difficult to describe the 
folding mechanisms in an unambiguous way. With the help of the atomic force microscopy, 
mechanical force has been used to prepare well defined initial states of proteins [H, [1]. Using 
the initial force, //, which is higher than the equilibrium critical force, fc, to unfold the 
tandem of poly ubiquitin (Ub), Fernandez and Li j2| have shown that the refolding can be 
initiated starting from stretched conformations or force denaturated ensemble (FDE) and 
quenching the force to a low constant value, fq {fq < fc)- Monitoring folding events as a 
function of the end-to-end distance (R) they have made the following important observations. 
1) Contrary to the standard folding from the thermal denaturated ensemble (TDE) the 
refolding under the quenched force is a multiple stepwise process and 2) The force-quench 
refolding time obeys the Bell formula [3|, Tp ~ Tp exp{fqAx p / ksT) , where Tp is the folding 
time in the absence of the quench force and Axp is the average location of the transition 
state (TS). 

Motivated by the experiments of Ferandez and Li [2|, we have studied 0] the refolding 
of the domain 127 of the human muscle protein using the Cq,-Go model |5| and the four- 
strand /5-barrel model sequence SI [6| (for this sequence the nonnative interactions are also 
taken into account). Basically, we have reproduced qualitatively the major experimental 
findings listed above. In addition we have shown that the refolding is two-state process in 
which the folding to the native basin attractor (NBA) follows the quick collapse from initial 
stretched conformations with low entropy. The corresponding kinetics can be described 
by the bi-exponential time dependence, contrary to the single exponential behavior of the 
folding from the TDE with high entropy. 

In order to make the direct comparison with the experiments of Fernandez and Li j^, 
in this paper we performed simulations for a single domain Ub using the Cq-Go model (see 
Material and Method for more details). Because the study of refolding of 76-residue Ub (Fig. 
[T]a) by all-atom simulations is beyond present computational facilities the Go modeling is 
an appropriate choice. Most of the simulations have been carried out at T = 0.85Tp = 285 
K. Our present results for refolding upon the force quench are in the qualitative agreement 
with the experimental findings of Fernandez and Li, and with those obtained for 127 and 
SI theoretically A number of quantitative differences between 127 and Ub will be also 
discussed. For Ub we have found the average location of the transition state Axp ~ 0.96 
nm which is in reasonable agreement with the experimental value 0.8 nm 0]. 

Experimentally, the unfolding of the polyubiquitin has been studied by applying a con- 
stant force 0]. The mechanical unfolding of Ub has previously investigated using Go-like 
and all-atom models jl, 0]. In particular, Irback et al. have explored mechanical unfolding 
pathways of structures A, B, C, D and E (see the definition of these structures and the 
/3-strands in the caption to Fig. [T]) and the existence of intermediates in detail. We present 
our results on mechanical unfolding of Ub for five following reasons. First, the barrier to the 
mechanical unfolding has not been computed. Second, experiments of Schlierf et al. [3| have 
suggested that cluster 1 (strands SI, S2 and the helix A) unfolds after cluster 2 (strands 
S3, S4 and S5). However, this observation has not yet been studied theoretically. Third, 
since the structure C, which consists of the strands SI and S5, unzips first, Irback et al. 
pointed out that the strand S5 unfolds before S2 or the terminal strands follows the unfold- 
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ing pathway SI ^ S5 — S2. This conclusion may be incorrect because it has been obtained 
from the breaking of the contacts within the structure C. Fourth, in pulhng and force-clamp 
experiments the external force is applied to one end of proteins whereas the other end is 
kept fixed. Therefore, one important question emerges is how fixing one terminus affects the 
unfolding sequencing of Ub. This issue has not been addressed by Irback et al. Fifth, 
using a simplified all-atom model it was shown that mechanical intermediates occur more 
frequently than in experiments It is relevant to ask if a Cq,-Go model can capture similar 
intermediates as this may shed light on the role of non-native interactions. 

In this paper, from the force dependence of mechanical unfolding times we estimated the 
distance between the native state and the transition state to be ^xuf ~ 0.24 nm which 
is close to the experimental results of Carrion- Vazquez et al. 10| and Schlierf et al. 0]. 



In agreement with the experiments 0], cluster 1 was found to unfold after cluster 2 in our 
simulations. Applying the force to the both termini, we studied the mechanical unfolding 
pathways of the terminal strands in detail and obtained the sequencing SI ^ S2 ^ S5 which 
is different from the result of Irback et al.. When the N-terminus is fixed and the C-terminus 
is pulled by a constant force the unfolding sequencing was found to be very different from 
the previous case. The unzipping initiates, for example, from the C-terminus but not from 
the N-one. Anchoring the C-end is shown to have a little effect on unfolding pathways. 
We have demonstrated that the present Cq-Go model does not capture rare mechanical 
intermediates, presumably due to the lack of non-native interactions. Nevertheless, it can 
correctly describe the two-state unfolding of Ub 0]. 

It is well known that thermal unfolding pathways may be very different from the mechan- 
ical ones, as has been shown for the domain 127 This is because the force is applied 
locally to the termini while thermal fluctuations have the global effect on the entire pro- 
tein. In the force case unzipping should propagate from the termini whereas under thermal 
fluctuations the most unstable part of a polypeptide chain unfolds first. 

The unfolding of Ub under thermal fluctuations was investigated experimentally by 
Cordier and Grzesiek 1^ and by Chung et al. [l^]. If one assumes that unfolding is the 
reverse of the refolding process then one can infer information about the unfolding pathways 
from the experimentally determined 0- values 14] and values [l^, 1^. The most com- 
prehensive </)-value analysis is that of Went and Jacskon. They found that the C-terminal 
region which has very low 0-values unfolds first and then the strand SI breaks before full 
unfolding of the a helix fragment A occurs. However, the detailed unfolding sequencing of 
the other strands remains unknown. 

Theoretically, the thermal unfolding of Ub at high temperatures has been studied by 
all-atom molecular dynamics (MD) simulations by Alonso and Daggett 17j and Larios et 
al. [3]. In the latter work the unfolding pathways were not explored. Alonso and Daggett 
have found that the a-helix fragment A is the most resilient towards temperature but the 
structure B breaks as early as the structure C. The fact that B unfolds early contradicts not 
only the results for the 0-values obtained experimentally by Went and Jackson IJ] but also 
findings from a high resolution NMR 1^. Moreover, the sequencing of unfolding events for 
the structures D and E was not studied. 

What information about the thermal unfolding pathways of Ub can be inferred from 
the folding simulations of various coarse-grained models? Using a semi-empirical approach 
Fernandez predicted [l9| that the nucleation site involves the /5-strands SI and S5. This 
suggests that thermal fluctuations break these strands last but what happens to the other 
parts of the protein remain unknown. Furthermore, the late breaking of S5 contradicts the 
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unfolding l3] and folding 14 1 experiments. From later folding simulations of Fernandez et 



al. [20, l21| one can infer that the structures A, B and C unzip late. Since this information 
is gained from 0-values, it is difficult to determine the sequen cing of unfolding events even 
for these fragments. Using the results of Gilis and Rooman [22| we can only expect that 
the structures A and B unfold last. In addition, with the help of a three-bead model it 
was found [23| that the C-terminal loop structure is the last to fold in the folding process 
and most likely plays a spectator role in the folding kinetics. This implies that the strands 
S4, S5 and the second helix (residues 38-40) would unzip first but again the full unfolding 
sequencing can not be inferred from this study. 

Thus, neither the direct MD [13] nor indirect folding simulations 19, 20, 21, [22, 23l 
provide a complete picture of the thermal unfolding pathways for Ub. One of our aims is 
to decipher the complete thermal unfolding sequencing and compare it with the mechanical 
one. The mechanical and thermal routes to the denaturated states have been found to be 
very different from each other. Under the force the /3-strand SI, e.g., unfolds first, while 
thermal fluctuations detach strand S5 flrst. The later observation is in good agreement with 
NMR data of Cordier and Grzesiek 
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A detailed comparison with available experimental 
and simulation data on the unfolding sequencing will be presented. The free energy barrier 
to thermal unfolding was also calculated. 

To summarize, in this paper we have obtained the following novel results. We have shown 
that the refolding of Ub is a two-stage process in which the "burst" phase exists on very 
short time scales. The construction of the T — f phase diagram allows us to determine the 
equilibrium critical force fc separating the folded and unfolded regions. Using the exponen- 
tial dependence of the refolding and unfolding times on /, ^xp and ^xjjf were computed. 
Our results for /c, Axp and ^xjjp are in acceptable agreement with the experiments. It has 
been demonstrated that flxing the N-terminus of Ub has much stronger effect on mechanical 
unfolding pathways compared to the case when the C-end is anchored. In comparison with 
previous studies, we provide a more complete picture for thermal unfolding pathways which 
are very different from the mechanical ones. 



Materials and Methods 

Cq-Go model for Ub 

We use coarse-grained continuum representation for Ub in which only the positions of Cq,- 
carbons are retained. The interactions between residues are assumed to be Go-like and the 
energy of such a model is as follows [sl 
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Here Ac^j = 4>i — 4>oi, i^i,i+i is the distance between beads i and i + 1, 6i is the bond angle 
between bonds (i — 1) and i, and 0j is the dihedral angle around the ith bond and rij is the 
distance between the ith. and jth residues. Subscripts "0" , "NC" and "NNC" refer to the 
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native conformation, native contacts and non-native contacts, respectively. Residues i and 
j are in native contact if roij is less than a cutoff distance dc taken to be dc = 6.5 A, where 
roij is the distance between the residues in the native conformation. With this choice of dc 
and the native conformation from the PDB (Fig. [T]a), we have the total number of native 

contacts Qmax = 99. 

The first harmonic term in Eq. ([T]) accounts for chain connectivity and the second 
term represents the bond angle potential. The potential for the dihedral angle degrees of 
freedom is given by the third term in Eq. ([1]). The interaction energy between residues that 
are separated by at least 3 beads is given by 10-12 Lennard- Jones potential. A soft sphere 
repulsive potential (the fourth term in Eq. [1]) disfavors the formation of non-native contacts. 
The last term accounts for the force applied to C and N termini along the end-to-end vector 
R. We choose Kr = lOOe^/A^, Ke = 20eH/rad^, K^l^ = en, and K^f = O.ben, where en 
is the characteristic hydrogen bond energy and C = 4 A. Since Tp = O.Q75eH (see below) 
and Tp = 332.5K j2J], we have en = 4.1 kJ/mol = 0.98 kcal/mol. Then the force unit 
[/] = en/ A = 68.0 pN. 

We assume the dynamics of the polypeptide chain obeys the Langevin equation. The 



equations of motion (see Ref. l25| for details) were integrated using the velocity form of the 



Verlet algorithm [26| with the time step At = O.OOSr^, where Tp = {ma?' / euY^'^ ~ 3 ps. 
Simulations 

In order to obtain the T — f phase diagram we use the fraction of native contacts or the 



overlap function [27 

X 
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where Aij is equal to 1 if residues i and j form a native contact and otherwise, and 6{x) 
is the Heaviside function. The argument of this function guarantees that a native contact 
between i and j is classified as formed when Vij is shorter than 1.2roij. The probability of 
being in the native state, /at, which can be measured by various experimental techniques, is 
defined as Jn =< X >5 where < ... > stands for a thermal average. The T — f phase diagram 
( a plot of 1 — /at as a function off and T) and thermodynamic quantities were obtained by 



the multiple hi stog ram method [28| extended to the case when the external force is applied 
to the termini |29l . l30|. In this case the reweighting is carried out not only for temperature 
but also for force. We collected data for six values of T at / = and for five values of / 
at a fixed value of T. The duration of MD runs for each trajectory was chosen to be long 
enough to get the system fully equilibrated (9xl0^rL from which LSxlO^r^, were spent on 
equilibration). For a given value of T and / we have generated 40 independent trajectories 
for thermal averaging. 

For the mechanical unfolding we have considered two cases. In the first case the external 
force is applied via both termini N and C. In the second case it is applied to either N- or 
C-terminus. 

To simulate the mechanical unfolding the computation has been performed at T = 285 K 
and mainly at the constant force / = 70, 100, 140 and 200 pN. This allows us to compare our 
results with the mechanical unfolding experiments 0] and to see if the unfolding pathways 
change at low forces. Starting from the native conformation but with different random 
number seeds the unfolding sequencing of helix A and five /^-stands is studied by monitoring 
fraction of native contacts as a function of the end-to-end extension. In the case of structures 
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A, B, C, D and E we consider not only the evolution of the number of intra-structure contacts 
as has been done by Irback et al. [9(], but also the evolution of all contacts (intra-structure 
contacts and the contacts formed by a given structure with the rest of a protein). 

In the thermal unfolding case the simulation is also started from the native conformations 
and it is terminated when all of the native contacts are broken. Due to thermal fluctuations 
there is no one-to-one correspondence between R and time. Therefore R ceases to be a 
good reaction coordinate for describing unfolding sequencing. To rescue this, for each i-th 
trajectory we introduce the progressive variable 5i = t/r^p, where r^p is the unfolding time. 
Then we can average the fraction of native contacts over a unique window < 5^ < 1 and 
monitor the unfolding sequencing with the help of the progressive variable 5. 



Results 

Temperature-force phase diagram and thermodynamic quantities 

The T — f phase diagram, obtained by the extended histogram method (see Materials and 
Methods), is shown in Fig. [2]a. The folding-unfolding transition, defined by the yellow region, 
is sharp in the low temperature region but it becomes less cooperative (the fuzzy transition 
region is wider) as T increases. The weak reentrancy (the critical force slightly increases 
with T) occurs at low temperatures. This seemingly strange phenomenon occurs as a result 
of competition between the energy gain and the entropy loss upon stretching. The similar 
cold unzipping transition was also observed in a number of models for heteropolymers [sH 
and proteins [29| including the Co-Go model for 127 (MS Li, unpublished results). As follows 
from the phase diagram, at T = 285 K the critical force fc ~ 30 pN which is close to fc ~ 25 
pN, estimated from the experimental pulling data (To estimate fc from experimental pulling 
data we use fmax ~ /cln(w/fmjn) [12], where fmax is the maximal force needed to unfold a 
protein at the pulling speed v. From the raw data in Fig. 3b of Ref. [lo| we obtain ~ 25 
pN). Given the simplicity of the model this agreement can be considered satisfactory and it 
validates the use of the Go model. 

Figure [2^ shows the temperature dependence of population of the native state f^. Fitting 
to the standard two-state curve = , . „ \, — t ^ ,, ^i one can see that it works pretty 

•'^ l+exp[-A_ff,„(l-^)/fcsT] ^ 

well (solid curve) around the transition temperature but it gets worse at high T due to 
slow decay of /at. Such a behavior is characteristic for almost all of theoretical models 
2^ including the all-atom ones 33 1. In fitting we have chosen the hydrogen bond energy 



en = 0.98 kcal/mol in Hamiltonian ([T]) so that Tp = Tm = 0.675eH/kB coincides with the 
experimental value 332.5 K [13]. From the fit we obtain AiJm = 11.4 kcal/mol which is 
smaller than the experimental value 48.96 kcal/mol indicating that the Go model is, as 
expected, less stable compared to the real Ub. Taking into account non-native contacts and 
more realistic interactions between side chain atoms is expected to increase the stability of 
the system. 

The cooperativity of the denaturation transition may be characterized by the cooperativ- 



ity index, fi^ (see Refs. |3J and|35| for definition). From simulation data for /jv presented in 
Fig. [2f) we have Qc ~ 57 which is considerably lower than the experimental value Qc ~ 384 
obtained with the help of AH^ = 48.96 kcal/mol and Tm = 332. 5K [2^ . The underestima- 
tion of f2c in our simulations is not only a shortcoming of the off-lattice Go model |36| but 
also a common problem of much more sophisticated force fields in all-atom models |33| . 

Another measure of the cooperativity is the ratio between the van't Hoff and the calori- 
metric enthalpy K2 jl^]. For the Go Ub we obtained k,2 ~ 0.19. Applying the base line 
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subtraction 38| gives H2 ~ 0.42 which is still much below K2 ~ 1 for the truly one-or-none 
transition. Since ^2 is an extensive parameter, its low value is due to the shortcomings of 
the off-lattice Go models but not due to the finite size effects. More rigid lattice models give 
better results for the calorimetric cooperativity (39| . 

Figure [3]a shows the free energy as a function of Q for several values of force at T = Tp. 
Since there are only two minima, our results support the two-state picture of Ub 0, [llj . As 
expected, the external force increases the folding barrier, AFp (AFp = F^s — Fa) and it 
lowers the unfolding barrier, AFup {AFjjf = Fts — F^). From the linear fits in Fig. [3]& we 
obtain the average distance between the TS and D states, Axp = AFp/f ~ 1 nm, and the 
distance between TS and the native state, Axuf = AFup/f ~ 0.13 nm. Note that Axp 
is very close to Axf ~ 0.96 nm obtained from refolding times at a bit lower temperature 
T = 285 K (see Fig. [6] below). However, Axuf is lower than value 0.24 nm followed from 
mechanical unfolding data a^t f > (Fig. [H]). This difference may be caused by either 
sensitivity of Axjjp to the temperature or the determination of Axuf from the approximate 
free energy landscape as a function of a single coordinate Q is not sufficiently accurate. 

We have also studied the free energy landscape using i? as a reaction coordinate. The 
dependence of F on i? was found to be smoother (results not shown) compared to what was 
obtained by Kirmizialtin et al. using a more elaborated model [23] which involves the 



non-native interactions. 

Refolding under quenched force 

Our protocol for studying the refolding of Ub is identical to what has been done on the 
experiments of Fernandez and Li . We first apply the force // ~ 70 pN to prepare initial 
conformations (the protein is stretched if i? > 0.8L, where the contour length L = 28.7 nm). 
Starting from the FDE we quenched the force to fq < fc and then monitored the refolding 
process by following the time dependence of the number of native contacts Q(t), R(t) and 
the radius of gyration Rg{t) for typically 50 independent trajectories. 

Figure H] shows considerable diversity of refolding pathways. In accord with experiments 
[2I and simulations for 127 [3] , the reduction of R occurs in a stepwise manner. In the fq = 
case (Fig. |l]a) R decreases continuously from ^ 18 nm to 7.5 nm (stage 1) and fluctuates 
around this value for about 3 ns (stage 2). The further reduction to i? ^ 4.5 nm (stage 
3) until a transition to the NBA. The stepwise nature of variation of Q(t) is also clearly 
shown up but it is more masked for Rg{t). Although we can interpret another trajectory for 
fq = (Fig. Hb) in the same way, the time scales are different. Thus, the refolding routes 
are highly heterogeneous. 

The pathway diversity is also evident for fq>0 (Fig. Hie and d). Although the picture 
remains qualitatively the same as in the fq = case, the time scales for different steps 
becomes much larger. The molecule fluctuates around i? ~ 7 nm, e.g., for ^ 60 ns (stage 
2 in Fig. |l]c) which is considerably longer than ^ 3 ns in Fig. |l]a. The variation of Rg{t) 
becomes more drastic compared to the fq = case. 

Figure [5] shows the time dependence of < R(t) >,< Q(t) > and < Rg(t) >, where 
< ... > stands for averaging over 50 trajectories. The left and right panels correspond to 
the long and short time windows, respectively. For the TDK case (Fig. [5]a and b) the single 
exponential fit works pretty well for < R{t) > for the whole time interval. A little departure 
from this behavior is seen for < Q(t) > and < Rg(t) > for t < 2 ns (Fig. [5]6). Contrary 
to the TDK case, even for fq = (Fig. [5]c and d) the difference between the single and 
bi-exponential fits is evident not only for < Q(t) > and < Rg(t) > but also for < R(t) >. 
The time scales, above which two fits become eventually identical, are slightly different for 



7 



three quantities (Fig. [5](i). The failure of the single exponential behavior becomes more and 
more evident with the increase of fq, as demonstrated in Figs. [5]e and / for the FDE case 
with fg = 6.25 pN. 

Thus, in agreement with our previous results, obtained for 127 and the sequence SI 
01, starting from FDE the refolding kinetics compiles of the fast and slow phase. The 
characteristic time scales for these phases may be obtained using a sum of two exponentials, < 
A{t) >= Ao + Ai exp{-t / rf) + A2exp{-t/T^), where A stands for R, Rg or Q. Here 
characterizes the burst-phase (first stage) while may be either the collapse time (for R 
and Rg) or the folding time (for Q) (rf < r^). As in the case of 127 and SI 0], and 
Ti ^ are almost independent on fq (results not shown). We attribute this to the fact that 
the quench force {f^°'^ ~ 9 pN) is much lower than the entropy force (/e) needed to stretch 
the protein. At T = 285 K, one has to apply fe ~ 140 pN for stretching Ub to 0.8 L. 
Since f^""^ << /e the initial compaction of the chain that is driven by fe is not sensitive 
to the small values of fg. Contrary to rf, was found to increase with fg exponentially. 
Moreover, < ^2 " < implying that the chain compaction occurs before the acquisition 
of the native state. 

Figure [6] shows the dependence of the folding times on fg. Using the Bell-type formula 
jsl and the linear fit in Fig. [6] we obtain Aa^ ^ 0.96 nm which is in acceptable agreement 
with the experimental value Axp ~ 0.8 nm [2[]. The linear growth of the free energy barrier 
to folding with fg is due to the stabilization of the random coil states under the force. Our 
estimate for Ub is higher than Axp ~ 0.6 nm obtained for 127 One of possible reasons 
for such a pronounced difference is that we used the cutoff distance dc = 0.65 and 0.6 nm in 
the Go model ([T]) for Ub and 127, respectively. The larger value of dc would make a protein 
more stable (more native contacts) and it may change the free energy landscape leading to 
enhancement of Ax p. This problem requires further investigation. 

Absence of mechanical unfolding intermediates in Ca-Go model 

In order to study the unfolding dynamics of Ub, Schlierf et al. [3] have performed the AFM 
experiments at a constant force / = 100, 140 and 200 pN. The unfolding intermediates 
were recorded in about 5% of 800 events at different forces. The typical distance between 
the initial and intermediate states is AR = 8.1 ± 0.7 nm [3]. However, the intermediates 
do not affect the two-state behavior of the polypeptide chain. Using the all-atom models 
Irback et al. j9| have also observed the intermediates in the region 6.7 nm < R < 18.5 nm. 
Although the percentage of intermediates is higher than in the experiments, the two-state 
unfolding events remain dominating. To check the existence of force-induced intermediates 
in our model, we have performed the unfolding simulations for / = 70, 100, 140 and 200 pN. 
Because the results are qualitatively similar for all values of force, we present / = 100 pN 
case only. 

Figure [7] shows the time dependence of R{t) for fifteen runs starting from the native value 
Rn ~ 3.9 nm. For all trajectories the plateau occurs a.t R 4.4 nm. As seen below, passing 
this plateau corresponds to breaking of intra-structure native contacts of structure C. At 
this stage the chain ends get almost stretched out, but the rest of the polypeptide chain 
remains native-like. The plateau is washed out when we average over many trajectories 
and < R{t) > is well fitted by a single exponential (Fig. [7]), in accord with the two-state 
behavior of Ub Q • 

The existence of the plateau observed for individual unfolding events in Fig. [7| agrees 
with the all-atom simulation results of Irback et al. 191 who have also recorded the similar 
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plateau at R ~ 4.6 nm at short time scales. However unfolding intermediates at larger 
extensions do not occur in our simulations. This is probably related to neglect of the non- 
native interactions in the Cq,-Go model. Nevertheless, this simple model provides the correct 
two-state unfolding picture of Ub in the statistical sense. 

Mechanical unfolding barrier 

We now try to determine the barrier to the mechanical unfolding from the dependence of 
the unfolding times tuf on /. It should be noted that this way of determination of the 
unfolding barrier is exact and it would give a more reliable estimate compared to the free 
energy landscape approach in which the free energy profile is approximated as a function of 
only one order parameter. 

We first consider the case when the force is applied via both termini N and C. Since 
the force lowers the unfolding barrier, Tup should decrease as / increases (Fig. [S]). The 
present Go model gives tuf smaller than the experimental values by about eight orders of 
magnitude. E.g., for / = 100 pN, tuf ~ 12 ns whereas the experiments gives tuf ~ 2.77 s 
As seen from Fig. [8|, for / < 140 pN tuf depends on / exponentially. In this regime 
Tuf ~ TuF^^Pif^up/kBT), where Axuf is the average distance between the N and TS 
states. From the linear fit in Fig. [8] we obtained Axuf ~ 0.24 nm. Using different fitting 
procedures Schlierf et al. 0] have obtained Axuf ~ 0.14 nm and 0.17 nm. The larger value 
AxuF ~ 0.25 nm was reported in the earlier experiments Thus, given experimental 



uncertainty, the Cq-Go model provides a reasonable estimate of Axuf for the two-state Ub. 

In the high force regime (/ > 140 pN) instead of the exponential dependence tuf scales 
with / linearly (inset in Fig. [HD. The crossover from the exponential to the linear behavior is 
in full agreement with the earlier theoretical prediction |32(]. The similar crossover has been 
also observed [4l| for the another Go-like model of Ub but Axuf has not been estimated. 
At very high forces tuf is expected to be asymptotically independent of /. 

One can show that fixing one terminus of a protein has the same effect on unfolding 
times no matter the N- or C-terminus is fixed. Therefore, we show the results obtained 
for the case when the N-end is anchored. As seen from Fig. [HI the unfolding process is 
slowed down nearly by a factor of 2. It may imply that diffusion-collision processes 4^ 



play an important role in the Ub unfolding. Namely, as follows from the diffusion-collision 
model, the time, required for formation (breaking) contacts, is inversely proportional to the 
diffusion coefficient, -D, of a pair of spherical units. If one of them is idle, D is halved and 
the time needed to break contacts increases accordingly. Although fixing one end increases 
the unfolding times, it does not change the distance between the TS and the native state, 
Axuf (Fig. Ej. 

Mechanical unfolding pathways: force is applied to both termini 

Here we focus on the mechanical unfolding pathways by monitoring the number of native 
contacts as a function of the end-to-end extension AR = R — i?eq, where i?eq is the equi- 
librium value of R. For T = 285 K Rcq ~ 3.4 nm. Following Schlierf et al. 0], we first 
divide Ub into two clusters. Cluster 1 consists of strands SI, S2 and the helix A (42 native 
contacts) and cluster 2 - strands S3, S4 and S5 (35 native contacts). The dependence of 
fraction of intra-cluster native contacts is shown in Fig. [9] for / = 70 and 200 pN (similar 
results for / = 100 and 140 pN are not shown). In agreement with the experiments 0] the 
cluster 2 unfolds first. The unfolding of these clusters becomes more and more synchronous 
upon decreasing /. At / = 70 pN the competition with thermal fluctuations becomes so 
important that two clusters may unzip almost simultaneously. Experiments at low forces 
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are needed to verify this observation. 

The arrow in Fig. [9] marks the position AR = 8.1 nm, where some intermediates were 
recorded in the experiments 0] • At this point there is intensive loss of native contacts of the 
cluster 2 suggesting that the intermediates observed on the experiments are conformations 
in which most of the contacts of this cluster are already broken but the cluster 1 remains 
relatively structured (~ 40% contacts). One can expect that the cluster 1 is more ordered 
in the intermediate conformations if the side chains and realistic interactions between amino 
acids are taken into account. 

To compare the mechanical unfolding pathways of Ub with the all-atom simulation results 
jol we discuss the sequencing of helix A and structures B, C, D and E in more detail. We 
monitor the intra-structure native contacts and all contacts separately. The later include 
not only the contacts within a given structure but also the contacts between it and the rest 
of the protein. It should be noted that Irback et al. have studied the unfolding pathways 
based on the evolution of the intra-structure contacts. Fig. [TUk shows the dependence of the 
fraction of intra-structure contacts on AR at / = 100 pN. At AR f» Inm, which corresponds 
to the plateau in Fig. [3, most of the contacts of C are broken. In agreement with the all- 
atom simulations , the unzipping follows C— *B^D— >E— >A. Since C consists of the 
terminal strands SI and S5, it was suggested that these fragments unfold first. However, 
this scenario may be no longer valid if one considers not only intra-structure contacts but 
also other possible ones (Fig. [TOb). In this case the statistically preferred sequencing is B 
— >C— >-D^E— i>A which holds not only for /=100 pN but also for other values of /. If 
it is true then S2 unfold even before S5. To make this point more transparent, we plot the 
fraction of contacts for SI, S2 and S5 as a function of AR (Fig. [TTb) for a typical trajectory. 
Clearly, S5 detaches from the core part of a protein after S2 (see also the snapshot in Fig. 
[TTb). So, instead of the sequencing SI — > S5 — > S2 proposed by Irback et al, we obtain SI 
^ S2 ^ S5. 

The dependence of the fraction of native contacts on AR for individual strands is shown 
in Fig. [Hx (/ = 70 pN) and Fig. ^ (/=200 pN). At A = 8.1 nm contacts of SI, S2 and 
S5 are already broken whereas S4 and A remain largely structured. In terms of /3-strands 
and A we can interpret the intermediates observed in the experiments of Schlierf et al. 0] 
as conformations with well structured S4 and A, and low ordering of S3. This interpretation 
is more precise compared to the above argument based on unfolding of two clusters because 
if one considers the average number of native contacts, then the cluster 2 is unstructured in 
the intermediate state (Fig. [9]), but its strand S4 remains highly structured (Fig. [T2l) . 

From Fig. [12] we obtain the following mechanical unfolding sequencing 

SI ^ S2 ^ S5 ^ S3 ^ S4 ^ A. (3) 

It should be noted that the sequencing ([3]) is valid in the statistical sense. In some trajectories 
S5 unfolds even before SI and S2 or the native contacts of SI, S2 and S5 may be broken at 
the same time scale (Table 1). From the Table 1 it follows that the probability of having 
SI unfolded first decreases with lowering / but the main trend ([3]) remains unchanged. One 
has to stress again that the sequencing of the terminal strands SI, S2 and S5 given by Eq. [3] 
is different from what proposed by Irback et al. based on the breaking of the intra-structure 
contacts of C. Unfortunately, there are no experimental data available for comparison with 
our theoretical prediction. 

Mechanical unfolding pathways: One end is fixed 
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N-terminus is fixed. Here we adopted the same procedure as in the previous section except 
the N-terminus is held fixed during simulations. As in the process where both of the termini 
are subjected to force, one can show that the cluster 1 unfolds after the cluster 2 (results 
not shown). 

From Fig. [T3] we obtain the following unfolding pathways 

C^D^E^B^A, (4a) 

S5 ^ S3 ^ S4 ^ SI ^ S2 ^ A, (4b) 

which are also valid for the other values of force (/=70, 100 and 140 pN). Similar to the case 
when the force is applied to both ends, the structure C unravels first and the helix A remains 
the most stable. However, the sequencing of B, D and E changes markedly compared to the 
result obtained by Irback et al [9| (Fig. [TOk). 

As evident from Eqs. [3landl4b| anchoring the first terminal has a much more pronounced 
effect on the unfolding pathways of individual strands. In particular, unzipping commences 
from the C-terminus instead of from the N-one. Fig. [T3b shows a typical snapshot where one 
can see clearly that S5 detaches first. At the first glance, this fact may seem trivial because 
S5 experiences the external force directly. However, our experience on unfolding pathways of 
the well studied domain 127 from the human cardiac titin, e.g., shows that it may be not the 
case. Namely, as follows from the pulling experiments |43| and simulations El, the strand 



A from the N-terminus unravels first although this terminus is kept fixed. From this point 
of view, what strand of Ub detaches first is not a priori clear. In our opinion, it depends on 
the interplay between the native topology and the speed of tension propagation. The later 
factor probably plays a more important role for Ub while the opposite situation happens 
with 127. One of possible reasons is related to the high stability of the helix A which does 
not allow either for the N-terminal to unravel first or for seriality in unfolding starting from 
the C-end. 

C-terminus is fixed. One can show that unfolding pathways of structures A,B, C, D and 
E remain exactly the same as in the case when Ub has been pulled from both termini (see 
Fig. [T0|) . Concerning the individual strands, a slight difference is observed for S5 (compare 
Fig. [T3H and Fig. [12]). Most of the native contacts of this domain break before S3 and S4, 
except the long tail at extension Ai? ^11 nm due to high mechanical stability of only one 
contact between residues 61 and 65 (the highest resistance of this pair is probably due to 
the fact that among 25 possible contacts of S5 it has the shortest distance |61 — 65| = 4 in 
sequence). This scenario holds in about 90% of trajectories whereas S5 unravels completely 
earlier than S3 and S4 in the remaining trajectories. Thus, anchoring C-terminus has much 
less effect on unfolding pathways compared to the case when the N-end is immobile. 

It is worth to note that, experimentally one has studied the effect of extension geometry on 
the mechanical stability of Ub fixing its C-terminus 1^ . The greatest mechanical strength 



(the longest unfolding time) occurs when the protein is extended between N- and C-termini. 



This result has been supported by Monte Carlo 10| as well as MD [Sj simulations. However 
the mechanical unfolding sequencing has not been studied yet. It would be interesting to 
check our results on the effect of fixing one end on Ub mechanical unfolding pathways by 
experiments. 

Thermal unfolding pathways 

In order to study the thermal unfolding we follow the protocol described in Materials and 
Methods. Two hundreds trajectories were generated starting from the native conformation 
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with different random seed numbers. The fractions of native contacts of hehx A and five 
/3-strands are averaged over all trajectories for the time window < 6 < 1. The unfolding 
routes are studied by monitoring these fractions as a function of 6. Above T ^ 500 K the 
strong thermal fluctuations (entropy driven regime) make all strands and helix A unfold 
almost simultaneously. Below this temperature the statistical preference for the unfolding 
sequencing is observed. We focus on T = 370 and 425 K. As in the case of the mechanical 
unfolding the cluster 2 unfolds before cluster 1 (results not shown). However, the main 
departure from the mechanical behavior is that the strong resistance to thermal fluctuations 
of the cluster 1 is mainly due to the stability of strand S2 but not of helix A (compare Fig. 
[Hlc and (i with Fig. [T2l) . The unfolding of cluster 2 before cluster 1 is qualitatively consistent 
with the experimental observation that the C-terminal fragment (residues 36-76) is largely 
unstructured while native-like structure persists in the N-terminal fragment (residues 1-35) 



45l . I46l . |47|. This is also consistent with the data from the folding simulations [23!] as well 
as with the experiments of Went and Jackson [l3] who have shown that the 0-values ^ in 
the C-terminal region. However, our finding is at odds with the high 0-values obtained for 



several residues in this region by all-atom simulations [48|] and by a semi-empirical approach 
19| . One possible reason for high 0- values in the C-terminal region is due to the force fields. 
For example, Marianayagam and Jackson have em ploy ed the GROMOS 96 force field [i^ 
within the software GROMACS software package [50|. It would be useful to check if the 
other force fields give the same result or not. 

The evolution of the fraction of intra-structure contacts of A, B, C, D and E is shown 
in Fig. [T3]a (T = 425 K) and b (T =370 K). Roughly we have the unfolding sequencing, 
given by Eq. EH which strongly differs from the mechanical one. The large stability of the 
a helix fragment A against thermal fluctuations is consistent with the all-atom unfolding 
simulations [13] and the experiments [3]- The N-terminal structure B unfolds even after the 
core part E and at T = 370 K its stability is comparable with helix A. The fact that B can 
withstand thermal fluctuations at high temperatures agrees with the experimental results 



of Went and Jackson [ij] and of Cordier and Grzesiek [1^ who used the notation Pi/ (32 
instead of B. This also agrees with the results of Gilis and Rooman [221] who used a coarse- 
grained model but disagrees with results from all-atom simulations [iTj. This disagreement 
is probably due to the fact that Alonso and Daggett studied only two short trajectories 
and B did not completely unfold [l7| . The early unzipping of the structure C (Eq. [Sa|) is 
consistent with the MD prediction [l7|. Thus our thermal unfolding sequencing (Eq. [Sal) is 
more complete compared to the all-atom simulation and it gives the reasonable agreement 
with the experiments. 

We now consider the thermal unstability of individual /3-strands and helix A. At T = 370 
K (Fig. [T^ ) the trend that S2 unfolds after S4 is more evident compared to the T = 425 K 



case (Fig. [Mb). Overall, the simple Go model leads to the sequencing given by Eq. 

(C,D)^E^B^A (5a) 

S5^S3^S1^A^(S4,S2). (5b) 

From Eq. [31 Hb] and [5b] it is obvious that the thermal unfolding pathways of individual 
strands markedly differ from the mechanical ones. This is not surprising because the force 
should unfold the termini first while under thermal fluctuations the most unstable part 
is expected to detach first. Interestingly, for the structures the thermal and mechanical 
pathways (compare Eq. [5a] and [la]) are almost identical except that the sequencing of C and 
D is less pronounced in the former case. This coincidence is probably accidental. 
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The fact that S5 unfolds first agrees with the high-resolution NMR data of Cordier and 



Grzesiek [12| who studied the temperature dependence of hydrogen bonds of Ub. However, 
using the ip-value analysis Krantz et al [l^ have found that S5 (B3 in their notation) breaks 
even after SI and S2. One of possible reasons is that, as pointed out by Fersht jsij, if 
there is any plasticity in the transition state which can accomodate the crosslink between 
the metal and bi-histidines, then -^/^-values would be significantly greater than zero even for 
an unstructured region, leading to an overestimation of structure in the transition state. In 
agreement with our results, the </>- value analysis [3] yields that S5 breaks before SI and A 
but it fails to determine whether S5 breaks before S3. By modeling the amide I vibrations 
Chung et al. 11] argued that SI and S2 are more stable than S3, S4 and S5. Eq. [5b] shows 



that the thermal stability of SI and S2 is indeed higher than S3 and S5 but S4 may be 
more stable than SI. The reason for only partial agreement between our results and those of 
Chung et al. remains unclear. It may be caused either by the simplicity of the Go model or 
by the model proposed in Ref . l3] . The relatively high stability of S4 (Eq. [5b|) is supported 



by the ^z^- value analysis [1| 
Thermal unfolding barrier 

Figure [T5l shows the temperature dependence of the unfolding time tuf which depends on 
the thermal unfolding barrier, AF^p, exponentially, tuf ~ TuF^^vi^^uF/^BT). From the 
linear fit in Fig. [T5]we obtain AF^p ^ 10.48e;i ~ 10.3 kcal/mol. It is interesting to note that 
AF^p is compatible with AHm ~ 11.4 kcal/mol obtained from the equilibrium data (Fig. 
12^). However, the latter is defined by an equilibrium constant (the free energy difference 
between native and denatured states) but not by the rate constant (see, for example, Ref. 



52|). 



Discusion 

We have studied the refolding of Ub following the same protocol as in the force-clamp 
experiments of Fernandez and Li [2|. Under the low quenched force the refolding is two- 
stage process characterized by two different time scales rf^ and r^, where rf^ « r^. This 
result further strengthens our previous prediction [3] that the nature of the folding starting 
from the FDE does not depend on details of models. The simple C^-Go model provides 
reasonable estimates for the equilibrium critical force fc as well as the averaged distance 
between the D and TS states, Axp, and the distance between the N and TS states, Axup. 
We have also obtained AHm from the two-state fit of the population of the NBA, /tv, and the 
thermal unfolding barrier AF^p. It would be interesting to measure AF^p experimentally 
and compare it with AHm- 

The shortcoming of the Go model we used is its failure to capture seldom unfolding 
intermediates observed in the experiments 0] as well as in the all-atom simulations j^. 
However it mimics the overall two-state behavior of Ub. Our simulations suggest that the 
non-native interactions, neglected in the Go model, may be the cause of mechanical 
unfolding intermediates. 

Due to thermal fluctuations the thermal unfolding pathways are not well defined as in 
the mechanical case. Nevertheless, at T < 500 K the statistical preference in the sequencing 
of unfolding events is evident. In accord with the experiments, the cluster 2 unfolds before 
cluster 1 in the mechanical as well as in the thermal cases. However, in terms of individual 
strands we predict that mechanical and thermal unfolding follows very different pathways 
(Eq. [3] and Eq. [5bl) . Mechanically strand SI is the most unstable whereas the thermal 
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fluctuations break contacts of S5 first. If we consider only breaking of intra-structure native 
contacts, then our mechanical sequencing agrees with the all-atom simulation results j^. 
It is probably not unexpected because mechanical unfolding pathways may depend largely 
on the topology of the native conformation and in some cases the Go-like models may give 
results comparable with experimental ones [8|. However, contrary to Irback et al. p], we 
predict that the terminal strands follow the mechanical unfolding sequencing: SI — * S2 — >■ 
S5. It would be very exciting to perform the AFM experiments to verify this prediction and 
the whole unfolding sequencing (Eq. [3]) 

We have considered the effect of fixing one end on unfolding kinetics and found that it 
delays the unfolding by nearly a factor of 2 regardless to what end is anchored. We argue 
that this general result rnay be understood, using the diffusion-collision model developed 
by Karplus and Weaver [4^. However, fixing one terminus does not affect the distance 
between the native state and TS. One of the most interesting results is that what terminus 
we keep fixed matters for the unfolding sequencing. Namely, anchoring the N-end changes 
it dramatically (see Eq. [3] and Eq. l4bl) whereas fixing the C-end has only a minor effect. 

As evident from Eqs. Ea] and [5b] and the detailed discussion in the Introduction, our 
ihermal unfolding sequencing is more complete compared with previous theoretical studies 




231 ]. We have obtained some agreement with the experimental data 



15l | on the instability of the structures and /5-strands. However the picture for 



thermal unfolding pathways is still incomplete. More experiments are needed to check our 
prediction given by Eqs. [5a] and [5b] 

We have also shown that refolding from FDE and folding from TDE have the same 
pathways which are not sensitive to the quenched force. The refolding/folding sequencing is 
the same as for the thermal unfolding (see Eqs. [5aland [5bl) but in the inverse order implying 
that the protein folding is the reversible process. 
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TABLE 1. Dependence of unfolding pathways on the external force. There are three possi- 
ble scenarios: SI — S2 ^ S5, S5 ^ SI ^ S2, and three strands unzip almost simultaneously 
(S1,S2,S5). The probabilities of observing these events are given in percentage. 
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Figure Captions 



FIGURE 1. (a) Native state conformation of ubiquitin taken from the PDB (PDB ID: 
lubq). There are five /3-strands: SI (2-6), S2 (12-16), S3 (41-45), S4 (48-49) and S5 (65-71), 
and one hehx A (23-34). (b) Structures B, C, D and E consist of pairs of strands (S1,S2), 
(S1,S5), (S3,S5) and (S3,S4), respectively. In the text we also refer to helix A as the structure 
A. 

FIGURE 2. (a) The T — f phase diagram obtained by the extended histogram method. 
The force is applied to termini N and C. The color code for 1— < x(T, /) > is given on 
the right. The blue color corresponds to the state in the NBA, while the red color indicates 
the unfolded states. The vertical dashed line refers to T = 0.85Tp ^ 285 K at which most 
of simulations have been performed, (b) The temperature dependence of Jn (open circles) 
defined as the renormalized number of native contacts (see Material and Methods). The 
solid line refers to the two-state fit to the simulation data. The dashed line represents the 
experimental two-state curve with AHm = 48.96 kcal/mol and = 332. 5K |24| . 

FIGURE 3. (a) The dependence of the free energy on Q for selected values of / at T = Tp. 
D and N refer to the denaturated and native state, respectively, (b) The dependence of 
folding and unfolding barriers, obtained from the free energy profiles, on /. The linear fits 
y = 0.36 + 0.218x and y = 0.54 — 0.029a; correspond to AFp and AFjjp, respectively. From 
these fits we obtain Axp ~ 10 nm and Axuf ~ 0.13 nm. 

FIGURE 4. (a) and (b) The time dependence of Q, R and Rg for two typical trajectories 
starting from FDE (/g = and T = 285 K). The arrows 1, 2 and 3 in (a) correspond to 
time 3.1 {R = 10.9 nm), 9.3 {R = 7.9 nm) and 17.5 ns (i? = 5 nm). The arrow 4 marks the 
folding time Tp = 62 ns {R = 2.87 nm) when all of 99 native contacts are formed, (c) and 
(d) are the same as in (a) and (b) but for fq = 6.25 pN. The corresponding arrows refer to 
t = 7.5 {R = 11.2 nm), 32 {R = 9.4 nm), 95 ns {R = 4.8 nm) and Tp = 175 ns {R = 3.65 
nm). 

FIGURE 5. (a) The time dependence of < Q{t) >, < R{t) > and < Rg{t) > when the 
refolding starts from TDE. (b) The same as in (a) but for the short time scale, (c) and (d) 
The same as in (a) and (b) but for FDE with fg = 0. (e) and (f) The same as in (c) and 
(d) but for ^=6.25 pN. 

FIGURE 6. The dependence of folding times on the quench force at T = 285 K. Tp was 
computed as the average of the first passage times {tp is the same as extracted from the 
bi-exponential fit for < Q(t) >). The result is averaged over 30 - 50 trajectories depending 
on fq. From the linear fit y = 3.951 -|- 0.267x with correlation level equal -0.96 we obtain 
Xp ~ 0.96 nm. 

FIGURE 7. Time dependence of the end-to-end distance for / = 100 pN. The thin curves 
refer to fifteen representative trajectories. The averaged over 200 trajectory < R{t) > is 
represented by the thick line. The dashed curve is the single exponential fit < R{t) >= 
21.08 — 16.81 exp(—x/r(7i?), where tuf ~ 11-8 ns. 
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FIGURE 8. Dependence of mechanical unfolding time on the force. Circles refer to the 
process when the force is applied to both N and C termini. Squares signifies the case when 
the N-end is fixed and the C-end is pulled. For the first case the linear fit {y = 9.247— 0.067a;) 
gives the distance between the native state and TS Axuf ~ 0.24 nm. In the second case, 
from the linear fit {y = 9.510 — 0.062x) we obtained Axup ~ 0.22 nm. Thus, within error 
bars fixing one end does not affect the value of Axup. The inset shows the linear dependence 
of on / in the high force regime. 

FIGURE 9. The dependence of fraction of the native contacts on the end-to-end extension 
for cluster 1 (solid lines) and cluster 2 (dashed lines) at / = 70pN and 200 pN. The results 
are averaged over 200 independent trajectories. The arrow points to the extension AR = 
8.1 nm. 

FIGURE 10. (a) The dependence of fraction of the intra-structure native contacts on AR 
for structures A, B, C, D and E at / = 100 pN. (b) The same as in a) but for all native 
contacts. The results are averaged over 200 independent trajectories. 

FIGURE 11. (a) The dependence of fraction of the native contacts on AR for strand SI, 
S2 and S5 (/ = 200j9A^) . The vertical dashed line marks the position of the plateau at AR ^ 
1 nm. (b) The snapshot, chosen at the extension marked by the arrow in a), shows that S2 
unfolds before S5. At this point all native contacts of SI and S2 have already broken while 
50% of the native contacts of S5 are still present. 

FIGURE 12. (a) The dependence of fraction of the native contacts on extension for A and 
all /3-strands at / = 70pN. (b) The same as in (a) but for / = 200 pN. The arrow points to 
AR = 8.1 nm where the intermediates are recorded on the experiments 0]. The results are 
averaged over 200 trajectories. 

FIGURE 13. (a) The dependence of fraction of the intra-structure native contacts on 
extension for all structures at / = 200pA^. The N-terminus is fixed and the external force 
is applied via the C-terminus. (b) The same as in (a) but for the native contacts of all 
individual /9-strands and helix A . The results are averaged over 200 trajectories, (c) A 
typical snapshot which shows that S5 is fully detached from the core while Si and S2 still 
have ~ 50% and 100% contacts, respectively, (d) The same as in (b) but the C-end is 
anchored and N-end is pulled. The strong drop in the fraction of native contacts of S4 at 
AR ^7.5 nm does not correspond to the substantial change of structure as it has only 3 
native contacts in total. 

FIGURE 14. (a) The dependence of fraction of intra-structure native contacts on the 
progressive variable 6 for all structures at T=425 K. (b) The same as in (a) but for T = 370 
K. (c) The dependence of the all native contacts of the /3-strands and helix A at T=425 K. 
(d) The same as in (c) but for T = 370 K. 

FIGURE 15. Dependence of thermal unfolding time tuf on en/T, where en is the hydrogen 
bond energy. The straight line is a fit y = —8.01 + 10.48a;. 
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